class ASTroCAT extends Ecosystem {
  
  // the ASTroCAT ecosystem model itself.
  // 
  // major changes since v1.6 (the Ecol. Modelling 2011 paper):
  // no stochasticity support
  // detritus and light dependence made explicit
  
  Statevar Ptot, Ztot, N, SD, LD, PAR;
  Statevar[] P, Z;
  
  Flux[][] FgrazP; // one for each P-Z pair
  Flux[] Fegest, Fexcret; // summed over P for each Z
  Flux[] Fuptake, FcoagP, Fzmort;
  Flux FcoagSD, FreminSD, FreminLD;
  
  FreeParam PMinEsd, PMaxEsd, Pesd0, Zesd0, mu00, mu0exp, ks0, ksexp, I00, I0exp, xpreyopt0, xpreyoptexp; // allometric curves
  FreeParam alpha, attSW, attP, K0, dxprey, epsilon, fegest, Rcoag, zeta, Rremin; // other params
  float[] Pesd, Zesd, mu0, ks, I0, xpreyopt; // params that vary among P, Z classes 
  float[][] phiP; // grazing preferences P -> Z
  int NP, NZ;
  
  // system definition ----------------------
  
  ASTroCAT() {}
  
  ASTroCAT(int NP) {
    this.NP = NP;
    this.NZ = NP;
    
    PMinEsd = addParam("min phyto ESD", "PMinEsd", "µm", 1);
    PMaxEsd = addParam("max phyto ESD", "PMaxEsd", "µm", 20);
    
    // state variables
    N = addVar("nutrients","N","µM");
    P = new Statevar[NP];
    for (int i=0; i<NP; i++) P[i] = addVar("phytoplankton-" + i, "P-" + i, "µM");
    Z = new Statevar[NZ];
    for (int j=0; j<NZ; j++) Z[j] = addVar("zooplankton-" + j, "Z-" + j, "µM");
    SD = addVar("small detritus","SD","µM");
    LD = addVar("large detritus","LD","µM");
    Ptot = addGhostVar("total P","Ptot","µM");
    Ztot = addGhostVar("total Z","Ztot","µM");
    PAR = addGhostVar("PAR","PAR","W/m^2");
    
    // minimum biomass values: always protect a little bit of seed stock. I don't think this is working!
    for (int i=0; i<NP; i++) P[i].min = 0.01;
    for (int j=0; j<NZ; j++) Z[j].min = 0.01;
    
    // fluxes
    Fuptake = addFluxSet("uptake","uptake", N, P);
    FcoagP = addFluxSet("P coagulation","coagP", P, LD);
    FgrazP = addFluxSet2D("grazing on P","grazP", P, Z);
    Fegest = addFluxSet("egestion","egest", Z, SD);
    Fexcret = addFluxSet("excretion","excret", Z, N);
    Fzmort = addFluxSet("Z mortality","Zmort", Z, SD);
    FcoagSD = addFlux("SD coagulation","coagD", SD, LD);
    FreminSD = addFlux("SD remineralization","reminSD", SD, N);
    FreminLD = addFlux("LD remineralization","reminLD", LD, N);    

    // allometric params are defined using a scale size of 1 µm for easy reading; the scale size is adjusted later
    Pesd0 = addParam("P scale size", "Pesd0", "µm", 1, 1, 100);
    Zesd0 = addParam("Z scale size", "Zesd0", "µm", 1, 1, 100);

    // P growth params
    mu00 = addParam("P max growth","mu00","1/d", 2.6); // ... so far, at 1 µm. But this will be changed before the user sees this string
    mu0exp = addParam("allometric exponent on mu0", "mu0exp", "", -0.45);
    ks0 = addParam("nutrient half-saturation","ks0","µM", 0.1); // at 1 µm
    ksexp = addParam("allometric exponent on ks", "ksexp", "", 1);
    alpha = addParam("initial P-E slope","alpha","", 0.15, 0.01, 0.2);
    attSW = addParam("light attenuation by seawater","attSW","1/m",0.1, 0, 0.2); // the Environment is repsonsible for using these to calculate the light profile
    attP = addParam("light attenuation by phytoplankton","attP","1/µM/m", 0.02, 0, 0.05);
    
    // grazing params
    I00 = addParam( "Z max grazing","I00","1/d", 26); // at 1 µm
    I0exp = addParam("allometric exponent on I0", "I0exp", "", -0.4);
    K0 = addParam("grazing half-saturation","K0","µM", 3, 0.3, 30);
    K0.logScale = true;
    xpreyopt0 = addParam("optimal prey size", "xpreyopt0", "µm", 0.65); // at 1 µm
    xpreyoptexp = addParam("allometric exponent on xpreyopt", "xpreyoptexp", "", 0.56);
    dxprey = addParam("prey size preference width","dxprey","", 0.25, 0.01, 1);
    dxprey.logScale = true;
    addParam("max preference","maxPref","", 10, 0, 10);
    epsilon = addParam("Z growth efficiency","epsilon","", 0.33, 0, 1);
    fegest = addParam("egestion fraction","fegest","",0.33, 0, 1); // a fraction of total grazing
    
    // loss and regeneration params
    Rcoag = addParam("P, SD coagulation rate","Rcoag","1/µM/d", 0.1, 0.01, 1);
    Rcoag.logScale = true;
    zeta = addParam("Z mortality","zeta","1/µM/d", 1, 0, 4);
    zeta.highlight = true;
    Rremin = addParam("SD, LD remineralization rate","Rremin","1/µM/d", 0.1, 0.01, 1);
    Rremin.logScale = true;
    addParam("stress-based P coagulation on/off", "do_stressCoag", "", 0, 0, 1).quantize(1);

    // internal arrays for storing some useful things: populated by updatePrivateParams()
    Pesd = new float[NP];
    Zesd = new float[NZ];
    mu0 = new float[NP];
    ks = new float[NP];
    I0 = new float[NZ];
    xpreyopt = new float[NZ];
    phiP = new float[NP][NZ];

    // set a better scale size, so the monotype version doesn't run with 1 µm cells
    Pesd0.initial = sqrt(PMinEsd.initial*PMaxEsd.initial);
    Zesd0.initial = pow(Pesd0.initial / xpreyopt0.initial, 1.0 / xpreyoptexp.initial); // set Z sizes such that each P is optimal prey for one Z
    mu00.initial *= pow(Pesd0.initial, mu0exp.initial); // really, (Pesd0 / previous scale size)^mu0exp: assuming previous scale size is 1
    ks0.initial *= pow(Pesd0.initial, ksexp.initial);
    I00.initial *= pow(Zesd0.initial, I0exp.initial);
    xpreyopt0.initial *= pow(Zesd0.initial, xpreyoptexp.initial);
    // adjust the names of the allometric params so the user can see what the sliders control
    mu00.name += " at " + round(Pesd0.initial) + " µM";
    ks0.name += " at " + round(Pesd0.initial) + " µM";
    I00.name += " at " + round(Zesd0.initial) + " µM";
    xpreyopt0.name += " at " + round(Zesd0.initial) + " µM";   
   
    // only make selected params visible to the user
    for (int i=0; i<params.length; i++) params[i].hidden = true;
    Rcoag.hidden = false;
    Rremin.hidden = false;
    zeta.hidden = false;
    mu00.hidden = false;
    mu0exp.hidden = false;

  } // end of constructor
  
  
  // ------------------------------
  
  
  void updatePrivateParams() {
    // update phytoplankton size and size-dependent params
    for (int i=0; i<NP; i++) {
      Pesd[i] = exp(map(i, 0, NP-1, log(PMinEsd.current), log(PMaxEsd.current)));
      mu0[i] = mu00.current * pow(Pesd[i] / Pesd0.current, mu0exp.current);
      ks[i] = ks0.current * pow(Pesd[i] / Pesd0.current, ksexp.current);
    }
    // zooplankton size and size-dependent params
    for (int j=0; j<NZ; j++) {
      Zesd[j] = Zesd0.current * pow(Pesd[j] / xpreyopt0.current, 1.0/xpreyoptexp.current); // set Z sizes such that each P is optimal prey for one Z
      I0[j] = I00.current * pow(Zesd[j] / Zesd0.current, I0exp.current);
      xpreyopt[j] = xpreyopt0.current * pow(Zesd[j] / Zesd0.current, xpreyoptexp.current);
    }
    // grazing preferences
    for (int j=0; j<NZ; j++) {
      for (int i=0; i<NP; i++) {
        phiP[i][j] = exp(-sq((log10(Pesd[i])-log10(xpreyopt[j])) / dxprey.current));
      }
    }    
  }
  
  
  // ------------------------------

  float calcUptake(int i, int k) {
    float nutlim = N.current[k] / (ks[i] + N.current[k]);
    float lightlim = alpha.current * PAR.current[k] / sqrt(sq(mu0[i]) + sq(alpha.current * PAR.current[k]));
    return mu0[i] * nutlim * lightlim * P[i].current[k];
  }
  
  
  void calcFluxes() {
    for (int k=0; k<Ncells; k++) { // for each spatial cell...
      // phyto growth
      for (int i=0; i<NP; i++) {
        Fuptake[i].current[k] = calcUptake(i,k);
      }
      
      // grazing
      // --- assemble all diet elements
      float[][] PscaledAsFood = new float[NP][NZ];
      float[] totalFood = new float[NZ];
      for (int j=0; j<NZ; j++) {
        totalFood[j] = 0;
        for (int i=0; i<NP; i++) {
          PscaledAsFood[i][j] = phiP[i][j] / K0.current * P[i].current[k];
          totalFood[j] += PscaledAsFood[i][j];
        }
      }
      // --- initialize excretion and egestion fields
      for (int i=0; i<NP; i++) Fegest[i].current[k] = 0;
      for (int j=0; j<NP; j++) Fexcret[j].current[k] = 0;
      // --- calculate fluxes
      for (int j=0; j<NZ; j++) {
        for (int i=0; i<NP; i++) {
          FgrazP[i][j].current[k] = I0[j] * Z[j].current[k] * PscaledAsFood[i][j] / (1. + totalFood[j]);
          float feg = (1-epsilon.current) * Pesd[i] / (5. + Pesd[i]); // egestion fraction is a function of Pesd, larger for cells >> 5 µm
          Fegest[j].current[k] += feg * FgrazP[i][j].current[k];
          Fexcret[j].current[k] += (1.0 - feg - epsilon.current) * FgrazP[i][j].current[k];          
        }
      }
      
      // coagulation
      for (int i=0; i<NP; i++) {
        FcoagP[i].current[k] = Rcoag.current * P[i].current[k];
        FcoagP[i].current[k] *= (Pesd[i] / (5. + Pesd[i])); // coagulation is a function of Pesd, larger for cells >> 5 µm
      }
      FcoagSD.current[k] = Rcoag.current * SD.current[k];
      
      // Z mortality
      for (int j=0; j<NZ; j++) {
        Fzmort[j].current[k] = zeta.current * Z[j].current[k] * Ztot.current[k];
      }
      
      // remineralization
      FreminSD.current[k] = Rremin.current * SD.current[k];
      FreminLD.current[k] = Rremin.current * LD.current[k];
      
      // diagnostics
      Ptot.current[k] = 0;
      for (int i=0; i<NP; i++) Ptot.current[k] += P[i].current[k];
      Ztot.current[k] = 0;
      for (int j=0; j<NZ; j++) Ztot.current[k] += Z[j].current[k];
    }
  }
  
  
  void saveSpecialInfo(Archive archive) {
    if (archive != null) {
      archive.addNameValuePair("mu0", mu0);
      archive.addNameValuePair("I0", I0);
      archive.addNameValuePair("ks", ks);
      archive.addNameValuePair("xpreyopt", xpreyopt);
      archive.addNameValuePair("phiP", phiP); 
      archive.addNameValuePair("Pesd", Pesd);
      archive.addNameValuePair("Zesd", Zesd);
    }
  }


}





// ----------------------------------------------------------------------------------------------------


class ASTroCAT_NetworkView extends NetworkView {
  
  ASTroCAT eco;
  
  ASTroCAT_NetworkView(Ecosystem eco) {
    this.eco = (ASTroCAT)eco;
  }
  
  void define() {
    varScale = 10 * (40.0/eco.NP);
    fluxScale = 0.5 * (40.0/eco.NP);

    float R = 0.8;
    float xcen = 0.2;
    float ycen = 0.1;
    int margin = 5;
    float arcInterval = TWO_PI/5.;
    float arcSpan = TWO_PI/5. * 0.8;

    float[] px = new float[eco.NP];
    float[] py = new float[eco.NP];
    float[] zx = new float[eco.NZ];
    float[] zy = new float[eco.NZ];
      
    float theta0 = degrees(4*arcInterval - PI/2.);
    addDot(eco.N, xcen+R*cos(theta0), ycen+R*sin(theta0), styles.dkBlueColor);
    theta0 = 3*arcInterval - PI/2.;
    float[][] vertices = fluxArc(xcen+R*cos(theta0+arcSpan/2.), ycen+R*sin(theta0+arcSpan/2.), xcen+R*cos(theta0-arcSpan/2.), ycen+R*sin(theta0-arcSpan/2.), -degrees(arcSpan), eco.NP);
    for (int i=0; i<eco.NP; i++) {
      px[i] = vertices[i][0];
      py[i] = vertices[i][1];
      Dot d = addDot(eco.P[i], px[i], py[i], transparency(styles.medGreenColor,0.8));
      d.name = nicenum(round(eco.Pesd[i]));
    }
    theta0 = 2*arcInterval - PI/2.;
    vertices = fluxArc(xcen+R*cos(theta0+arcSpan/2.), ycen+R*sin(theta0+arcSpan/2.), xcen+R*cos(theta0-arcSpan/2.), ycen+R*sin(theta0-arcSpan/2.), -degrees(arcSpan), eco.NZ);
    for (int j=0; j<eco.NZ; j++) {
      zx[j] = vertices[j][0];
      zy[j] = vertices[j][1];
      Dot d = addDot(eco.Z[j], zx[j], zy[j], transparency(styles.dkBrownColor,0.8));
      d.name = nicenum(round(eco.Zesd[j]));
    }
    theta0 = arcInterval - PI/2.;
    addDot(eco.SD, xcen+R*cos(theta0), ycen+R*sin(theta0), styles.medBlueColor);
    theta0 = - PI/2.;
    addDot(eco.LD, xcen+R*cos(theta0), ycen+R*sin(theta0), styles.ltBlueColor);
    
    // arrows
/*    Arrow[] aa = addArrowSet(eco.Fuptake, styles.ltGreenColor, 0);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.FgrazP, styles.dkBrownColor, 0);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.FcoagP, styles.dkGrayColor, 0);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.Fegest, styles.medGrayColor, 0);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.Fexcret, styles.medGrayColor, 0);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.Fzmort, styles.dkGrayColor, 0);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    addArrow(eco.FcoagSD, styles.dkGrayColor, 0);
    addArrow(eco.FreminSD, styles.ltBlueColor, 0);
    addArrow(eco.FreminLD, styles.ltBlueColor, 0);
*/
  }
  
  void updateArrows() {
    for (int i=0; i<arrows.length; i++) {
      arrows[i].update(fluxes[i].current[currentCell]);
    }
    if (arrows.length < 400) { // only draw the arrows if there aren't too many
      for (int i=0; i<arrows.length; i++) arrows[i].draw();
    }
  }
  
  boolean update() {
    updateArrows();
    boolean clicked = updateDots();
    
    pushStyle();
    Axes ax = dots[0].ax;
    textSize(16);
    textAlign(CENTER);
    float PP = 0, export = 0;
    for (int i=0; i<eco.NP; i++) PP += eco.getFlux("uptake-" + i).current[0];
    PP *= 6.6*12*365/1000 * 20; // the 20 is H
    text("" + round(PP) + " g C m-2 y-1", ax.left()+0.5*ax.width(), ax.top()+0.5*ax.height());
    
    popStyle();
    
    return clicked;
  }

  
}



